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Abstract: We present an analysis of the Sommerfeld enhancement to dark matter annihi- 
lation in the presence of an excited state, where the interaction inducing the enhancement 
is purely off-diagonal, such as in models of exciting or inelastic dark matter. We derive 
a simple and accurate semi-analytic approximation for the s-wave enhancement, which is 
valid provided the mass splitting between the ground and excited states is not too large, 
and discuss the cutoff of the enhancement for large mass splittings. We reproduce previ- 
ously derived results in the appropriate limits, and demonstrate excellent agreement with 
numerical calculations of the enhancement. We show that the presence of an excited state 
leads to generically larger values of the Sommerfeld enhancement, larger resonances, and 
shifting of the resonances to lower mediator masses. Furthermore, in the presence of a mass 
splitting the enhancement is no longer a monotonic function of velocity: the enhancement 
where the kinetic energy is close to that required to excite the higher state can be up to 
twice as large as the enhancement at zero velocity. 
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1. Introduction 

In particle physics, perturbation theory is commonly employed to calculate annihilation 
and scattering cross sections, with higher-order terms in the perturbative expansion being 
neglected. Provided the theory is not strongly coupled, this is generally a good approxima- 
tion for relativistic particles, but at low velocities and in the presence of a long-range force 
(classically, when the potential energy due to the long-range force is comparable to the 
particles' kinetic energy), the perturbative approach breaks down. In the nonrelativistic 
limit, the question of how the long-range potential modifies the cross section for short-range 
interactions can be formulated as a scattering problem in quantum mechanics, with signif- 
icant modifications to the cross sections occuring when the particle wavefunctions are no 
longer well approximated by plane waves (so the Born expansion is not well-behaved). The 
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deformation of the wavefunctions due to a Coulomb potential was calculated by Sommer- 
feld in HI, yielding a ~ 1/u enhancement to the cross section for short-range interactions 
(where the long-range behavior due to the potential can be factorized from the relevant 
short-range behavior). 

Recent measurements of the positron and electron cosmic ray spectra by the PAMELA, 
ATIC, Fermi and H.E.S.S experiments have observed a rise in the positron fraction starting 
at £^ ~ 10 GeV and extending (at least) up to ~ 100 GeV [^, as well as a broad excess 
in the total e"*" + e~ spectrum extending from several hundred GeV to O(TeV) energies 
i, I, g. Observations by WMAP §, | also suggest an excess in microwave emission from 
the inner Galaxy, termed the "WMAP Haze" Q, which has been attributed to synchrotron 
radiation from a new population of high energy electrons (with energies of tens to hundreds 
of GeV) |ll|. 



There has been considerable interest in these observations as possible signatures of 
dark matter (DM) annihilation or decay. If DM annihilation is responsible for the ex- 
cesses, the annihilation cross section must exceed that required for a WIMP initially in 
thermal equilibrium to freeze out with the correct relic density, by 1-4 orders of magnitude 
depending on the dark matter mass. Furthermore, to avoid antiproton bounds from the 
PAMELA experiment |12] and to generate sufficiently hard e~^e~ spectra, the dark matter 
should annihilate largely into leptons, pions, kaons and other light states (e.g. |13, |l4|). 

If dark matter is self-interacting, either via exchange of Standard Model gauge bosons 
or due to some novel interaction, then the nonperturbative "Sommerfeld enhancement" 
must be taken into account when computing annihilation cross sections in the present- 
day Galactic halo. The importance of the Sommerfeld enhancement in the context of dark 
matter annihilation has been studied by several authors [|^, [l^, 17, |l8|. Due to the velocity 
dependence of the Sommerfeld effect, and the presence of large resonances due to near- 
threshold bound states at low velocity (see e.g. [|l9|), the annihilation cross section could 
potentially be greatly increased in the present-day Galactic halo (/3 ~ 10~^ — 10~^), without 
greatly perturbing the freezeout of annihilations at /? ~ 0.3 (however, the Sommerfeld 
enhancement can modify the relic density at some level, particularly in the presence of 
resonances, see e.g. |2^, 21 1). 

Furthermore, it was pointed out in ||2^, 23| that the presence of light (< GeV) force car- 
riers coupling to the dark matter would naturally lead to unusual dark matter annihilation 
signatures, with zero or small branching ratios into heavy hadrons and gauge bosons, and 
hard spectra of leptons, pions and other light states. The dominant annihilation channel 
in such models would be annihilation of the dark matter into the new force carriers, which 
would then decay into kinematically accessible SM states. In particular, provided the mass 
of the new force carrier was less than twice the proton mass, no excess antiprotons would 
be produced. A new GeV-scale force in the dark sector could therefore naturally gener- 
ate both a large annihilation cross section via Sommerfeld enhancement, and the observed 



e+e spectra, without violating antiproton constraints |24, 25 



As noted in [24|, if dark matter is charged under some new dark gauge group, then 
when the new gauge boson acquires an O(GeV) mass the states in the dark matter mul- 
tiplet can naturally be split from each other by a small (©(MeV)) amount (this happens 
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automatically, due to loops of the new dark gauge bosons, if the gauge group is non-Abelian 



1 24, but such mass splittings can also be naturally generated in Abelian models, e.g. 
p7[] ). Furthermore, if the dark matter is a Majorana fermion or real scalar, its couplings 
with any vector boson must be off-diagonal, since it cannot carry conserved charge. Models 
of this type, where the dark matter can only scatter inelastically at tree level (that is, by 
exciting some slightly more massive state) are also motivated by the iDM and XDM 

|22] scenarios, which explain, respectively, the anomalies observed by the DAMA/LIBRA 
m and INTEGRAL/SPI [H experiments. 

Despite the fairly generic presence of these excited states, most recent analyses of the 
Sommerfeld enhancement in the context of the cosmic-ray anomalies have treated only the 
case of degenerate dark matter states, under the presumption that the introduction of an 
excited state does not significantly modify the enhancement to annihilation, unless the mass 



splitting is large enough to cause the enhancement to cut off (see |24] for an argument to 
this effect) ^. While it is true that the general qualitative features of the enhancement are 
mostly unchanged unless the mass splitting is quite large, modifications to the resonance 
structure occur for smaller mass splittings, and should be taken into account in any detailed 
analysis of the enhancement. However, the addition of even a single mass splitting adds 
an extra parameter to the problem, thus making numerical studies significantly more time- 
consuming; the numerical calculation of the enhancement can also become badly unstable 
in some parts of parameter space. 

In §|2| we derive simple expressions to aid in computing the Sommerfeld enhancement in 
the general case where the Ith partial wave dominates annihilation, for an arbitrary poten- 
tial coupling N states. In §^ we discuss the general properties of the case where the dark 
matter has a single excited state and the interaction is purely off-diagonal, including the 
limit of zero mass splitting and a characterization of the parts of parameter space where the 
enhancement is negligible. In §^ we derive an approximate semi-analytic solution for the 
Schrodinger equation in the case of a two-state system with purely off-diagonal interaction, 
focusing on the case of s-wave annihilation. In §|| we present a simple analytic approx- 
imation for the Sommerfeld enhancement, strictly valid when the annihilation channels 
and matrix elements are the same for two particles in the excited state as in the ground 
state, and also applicable to the more general case (up to a simple rescaling) , provided the 
enhancement is large. We employ this approximation to discuss the new features present in 
the case with a non-zero mass splitting, and confirm our results by numerical simulations. 

2. The Sommerfeld enhancement for a general A^-state system 

The Sommerfeld enhancement results from the distortion of the two-body wavefunction 
away from a plane wave, in the presence of a long-range potential. Provided that the 
range of the interaction responsible for annihilation is much shorter than the range of 
the potential, we can approximate the annihilation interaction as a delta function at the 
origin (in relative coordinates). Then the Sommerfeld-enhanced annihilation rate can be 



^The related problem of coUisional excitation of these states, via scattering mediated by a long-range 
force, has been considered by several authors, e.g. P2, p2{. 
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determined by first solving the scattering problem (for the chosen initial plane wave), and 
then contracting the state vector at the origin with the matrix describing the annihilation 
rate for the various states, which we will denote T. The fraction of scatterings involving 
transitions to the various states can be determined similarly, by solving the scattering 
problem for an initial plane wave purely in one two-body state, and then reading off the 
amplitudes of the (purely outgoing) partial waves in the other states. 

In terms of Feynman diagrams, the Sommerfeld enhancement can be expressed as the 
result of resumming an infinite series of ladder diagrams. Several authors have shown that 
in the non-relativistic limit, this procedure is equivalent to solving the Schrodinger equation 
with an appropriate (in general, matrix- valued) potential |l6|, 33, 34 1. 




Let us assume a spherically symmetric, real-valued N x N matrix potential, and 
transform to relative coordinates. Then as usual, separation of variables yields the ra- 
dial Schrodinger equation for the two-particle wavefunction. 



Riir) = lil + l)Ri{r), (2.1) 



for the Ith partial wave, where Ri{r) is a regular real A^-component vector function, 
is the dark matter mass (twice the reduced mass), and we incorporate any mass splitting 
terms into the potential matrix V (we assume throughout that any mass splittings are much 
smaller than the dark matter mass). There are 2N linearly independent (LI) solutions to 
this equation in total, but the requirement that all components be regular at the origin 
fixes of the boundary conditions. Let us promote Ri{r) to an N x N matrix, such 
that its columns are LI regular vector solutions Ri{r). The functions Ri{r)ij have the 
asymptotics as r — )• cxd, 

{Ri{r))i.j ^ sin (^kir - Utt + {6i)ij^ {ni),j, (2.2) 

where {ni)ij and {Si)ij are real numbers, and ki = ^y (m^v)'^ — m^Vn{oo)/h, where we 
assume that V(oo) is diagonal (since the elements of V that do not vanish as r — t- cxd 
correspond to mass terms, and we can choose to work in the basis of mass eigenstates). 
We will write the general solution to the Schrodinger equation in the form ^'(r, 9) = 
^i^Q Pi (cos 9) Ri{r)Ai, where Ai is a complex A^-vector giving the coefficients of the real 
vector functions in Ri{r). 

The scattering solution to the multi-state Schrodinger equation can be written, 

^„(r, 9) = cne'^-' + /„(r, 9)e'^-' /r, (2.3) 

where the A^-vector c describes the initial incoming plane wave and / describes the outgoing 
scattered wave. Using the asymptotic expansion of e*'^^. 



^ikz _^ 1 y-/2Z + l)P,(cos6')(e''='^ - (-n'e-^'^n, (2.4) 
2ikr 
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for each I we can write 

= p^(cos^)^ ^(ei(fenr-ii^+(50„,) _ e-<^"'-i'^+('^'W))(n/)„j(A)j. (2.5) 

Comparing coefficients yields, 

{21 + l)c„/A;n = i-'^e~^(^')"^(nz)„,(A)i. (2.6) 



This step corresponds to imposing the conditions, from Eq. 2.3, that (1) the scattered 
wave is purely outgoing, and (2) the normalization of the nth component of the initial 
plane wave is set by c„. 

Let us define an x matrix {Mi)ij = e~^^^^^*i {ni)ij, describing the asymptotics of 
the basis solutions for each I, and an x matrix C^j = {cj)n/kn, where 61,62, ...,c^ 
are A^-vectors that form a basis for possible incoming states (in the simplest example, we 
could simply choose {cj)n = Sjn)- Let us also define an N x N matrix Ai such that the 
mth column of Ai is the coefficient vector Ai corresponding to the incoming state described 
by Cm- Then Eq. I.t can be written in matrix form as Ai = i^(2l + l)Mf^C. A basis of 



solutions to the Schrodinger equation is then given by the columns of the matrix, 

^ = Y^ Piicos e)i\2l + l)RiM-^C. (2.7) 



1=0 



As r — ?■ 0, Ri{r) ~ r', so the lowest-/ partial wave will dominate (at least to the extent 
that the short-range interaction of interest can be approximated as a contact interaction). 
Writing Ri{r) = xi{f)/fi the reduced Schrodinger equation has the standard form, 



— X: {r) = V(r) - m^v^ + xi{r), (2.8) 

and Ri{r — )• 0) = r'x^'^^Ho)/(^ + !)!• Consequently we can write. 



^s-..Ar = 0) = x'(0)Mo-iC, ^i{r ^ 0) = PiicosOy (^^^) X^'^^Km^r^C. 

(2.9) 

Multiplying the annihilation matrix by this matrix and its conjugate (^'tp^') win give the 
Sommerfeld-enhanced annihilation cross-section. Specifically, the mth diagonal element of 
the resulting N x N matrix will give the enhancement for the ingoing state described by 
the vector Cm (since that element corresponds to ^mT^m)- In the case of higher than 
s-wave modes, it no longer suffices to treat the annihilation interaction as point-like for 
determining the annihilation cross section, but the enhancement relative to the unperturbed 
annihilation cross section can be calculated in this simple framework by taking the limit 
as r — )• 0. 
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Consider the unperturbed case with V(r) = V(oo), in which case V is purely diagonal. 
Then xi has an exact solution in terms of Bessel functions, 



{Xl{'r))ij = Vr {{Ci)ijJi+i/2{kir) + {C2)ijYi+i/2{kir)) . 
Imposing the condition of regularity at the origin sets C2 = 0, and the boundary condition 



of Eq. 2.2 sets (Ci)jj = ^J^Tki/2{nl)ij, as well as fixing 5i = 0. Then as r — )• the unper- 
turbed radial solution has the asymptotic form {xi{r = {nl)ij^/r^J■Kki/2{kir/2y^^/'^/T{l+ 
3/2) = (ni)ijA;^"*"V'+V(2/ + 1)!!, and Mi = ni, which in turn yields (^'/(r 0))ij = 
PiicosOy {^{21 + I) /{I + 1)!) ((/ + l)!/(2/ + 1)!!) k\{ej\. 

Now comparing the perturbed case to the unperturbed result, the Sommerfeld en- 
hancement where the lih. partial wave dominates, for the ingoing state described by Cm, 
becomes. 



(2/ + l)!! V {(x('+i)(0)M-iC)t . r • (x('+i)(0)M-iC)}, 
(/ + !)! 



(2.10) 



j j ( kj kj ) ^ ^i^ijiCm)j 

In the case where A^=l, soc— 7-1,(7— t-I/Zc and M; — )■ e"**^', this expression reduces to 



S 



(2/ + l)!!x('+^H0) 



(/ + l)!A;'^ 



(2.11) 



in agreement with |33, 34 1. 

Using this result requires us to know the phase shifts described by M;, since while the 
enhancement factor is real by definition, there are several different phase shifts involved 
and ^^^r^' contains cross terms (in contrast to the A'^ = 1 case where the phases cancel 
trivially, and so can be set to any arbitrary value). The phase shifts can be determined by 
solving the Schrodinger equation for an appropriate basis of regular solutions. 

An alternate form for the enhancement can be obtained in terms of irregular solu- 
tions to the (reduced) Schrodinger equation (Eq. |2.8|) , rather than the regular solution 
x(r). Consider the matrix pi{r) satisfying the reduced Schrodinger equation, such that 
{pi)ij{r) — )■ Tije^'^^^ as r — )• 00 (i.e. a purely outgoing spherical wave), for T some constant 



(complex) matrix to be determined, and {pi)ij{r — )■ 0) 



(this is the correct asymptotic 



form for the /th partial wave at small r, provided the potential grows more slowly than 1/r^ 
at small r, so the /(/ + l)/r^ term dominates). Note that if ki is imaginary, the "outgoing 
wave" boundary condition is replaced by requiring that component of the wavefunction to 
be exponentially decaying as r — )• 00. 

Now consider the quantity Wi = pfx'i ~ p'l^Xi- It is easy to check from the Schrodinger 
equation that Wi is independent of r so long as the potential matrix V is symmetric. But 
writing xii^ — ^ 0) = r^^^x''^^^\'^) / {I + !)!> we have. 



Wi{r = Q) 



21 + 1 
(7TT)! 



X 



(0), 



(2.12) 
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Wi{r oo)ij = ^ f Tmie*'''"''/;:^ cos f /c^r - ^/vr + j "-mj 



.p-*('5i)mj„ .1. 



m 



^ ^ {■^l)m,j km ■ (^-l^) 



m 



Taking C to be the diagonal matrix with elements 1/kn (i.e. the incoming states 
Cn are the natural basis states), becomes the diagonal matrix with elements kn, i.e. 
{C~^)ab = kaSab- Then we can write, 

Wi{r ^ oo)ij = i'Y,Tmi{Mi)aj{C-^)ma = i^T^ Mi),j , (2.14) 

m 

and by the r-independence of Wi we have, 

T^C-'M, = ^-^^±^X^^^'H0) x('+^)(0)Mr^C = z'^±^r^. (2.15) 

The enhancement to annihilation for particles initially in the mth state, if the Ith partial 
wave dominates, can then be written in the form, 

Sm=i ^\i ' ^ ^ (2.16) 

V "'m / mm 

In the s-wave case and assuming T real and symmetric, we recover the usual expression for 
the Sommerfeld-enhanced annihilation cross-section in an A^-state system (see e.g. [[T^]). 
In the absence of any potential, we recover 5 = 1 as required (note the identity r(/-|-l/2) = 
(2/ - l)!!V^/2'). 

3. The enhancement in a two-state system with off-diagonal interaction 

Consider the case of a U(l) vector interaction coupling two states with some small mass 
splitting 5, where the interacting particles cannot carry conserved charge and thus the 
long-range interaction between the |1) and |2) states is purely off-diagonal. Then the 4x4 
potential matrix coupling the 2-particle states (|11), |12), |21), |22)) is block diagonal (the 
1 11) and 1 22) states only couple to each other, and the same for the |12) and |21) states). 
We will be interested in the Sommerfeld enhancement for a pair of particles initially in the 
ground state, so we wish to solve the Schrodinger equation with the 2x2 matrix potential 
which couples the |11) and |22) states, given by |24|, 



—hca 



where the first row corresponds to the ground state |11) 
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Define the dimensionless parameters = {v/c)/a, es = Y^2<5/m^/a, = {m^/m^)/a. 
Then rescaling the r coordinate by am^c/h, we can rewrite the Schrodinger equation as, 



i{i+i) 



r{r) = ^ r ^(,). (3.2) 

( _ e2 _y(^^ \ 
More generally, consider the matrix M = " ui+i') o o some ar- 

bitrary scalar potential V{r). The eigenvalues and eigenvectors of M are given by. 



^l+{V{r)/{ei/2))^ 





(3.3) 

There is a transition in the behavior of the eigenvalues and eigenvectors at V{r) ~ £5/2. 
For V{r) e|/2, a 45° rotation approximately diagonalizes M, with the eigenvalues and 
eigenvectors taking the form, 

A±-^-^^ + |±m ^i-^iri- (3.4) 

For 1^(7') ^ £5/2, on the other hand, M is already approximately diagonal, with eigenvalues 
and eigenvectors: 

A+ (-v + ^& + r2 ~' ^+ I 1 h V'- 

In these regimes, because there is an r-independent diagonalization of the Schrodinger 
equation, if </>^(r') = A±(r)(/>-|-(r), then (t)±{r)'ilj± is an approximate solution to the matrix 
ODE ip"(r) = Mip(r). In the intermediate transition regime, however, where V{r) ~ £5/2, 
the eigenvectors are clearly not constant with respect to r, and attempting to diagonalize 
the differential equation gives rise to additional (non-diagonal) terms depending on the 
derivatives of the diagonalizing matrix. We can regard the behavior of the solution in 
this intermediate region as a nonadiabatic transition between the eigenstates, of Rosen- 
Zener-Demkov type (a general class of nonadiabatic transitions related to models first 
studied by Rosen and Zener [^] and Demkov [^] ) . We define the transition radius r-p by 
Virr) = ej/2. 

3.1 The limit of zero mass splitting 

In the limit as eg — t- 0, rj^ — t- 00, and the potential matrix is always exactly diagonalizable 
by a 45° rotation: the Sommerfeld enhancement then follows directly from the enhance- 
ments (or suppressions) for the attractive and repulsive scalar potentials. In particular, if 
all elements of the annihilation matrix T are identical, then the enhancement to s-wave 
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annihilation for particles originating in the ground state is given by. 
5*11 = 




(p+{oo) 


-M 


<A-(oo) - (1 


;.+ (oo) 


(p-{oo) - (j) 


+ (oo) 


4>+{oo) + q 


;._(oo) 




(/)+(oo) + 0-(oo) 


(j)-{oo) - 0+(oo) 


4>~{oo) - (/'+(oo) 


(j)+{oo) + 0_(oo) 



c« 



where (/>_ and 0+ are, respectively, the solutions to the scalar Schrodinger equation with 
attractive and repulsive potentials, with <j)-{0) = 4>+{0) = 1 (for the case of / > the same 
argument holds, with the appropriate rescaling of the enhancement (Eq. p.l6D and redefi- 
nition of (j)±{0)). So in this limit we precisely recover the usual result for the enhancement, 
and the repulsive component does not contribute. 

3.2 Regimes where the enhancement is negligible 

Before deriving a detailed expression for the enhancement, it is useful to consider the limits 
where the enhancement must approach unity, in terms of the dimensionless parameters e^,, 
and e^. 

The characteristic Bohr energy of the Coulomb potential is given by ~ a^m^c^. The 
presence of a Yukawa cutoff can only reduce this energy. If this scale is small compared 
to the kinetic energy, the potential is irrelevant and the enhancement must be 0(1): this 
condition is equivalent to my,v'^ > a^m^(? <^ €y>l. If the range of the Yukawa potential 
is much smaller than the Bohr radius, then again the potential cannot significantly affect 
the wavefunction. This criterion becomes l/m^ < l/arriy-^ i.e. ^ 1- These conditions 
are familiar from the usual 5 = case (see e.g. |]2^ ). 

We might suppose that if the energy of the mass splitting were significantly larger than 
the kinetic energy of the system, the higher state could not be excited and no interaction 
would occur. However, if the potential energy is sufficiently large, it can allow excitations 
at small r: the relevant criterion is therefore whether the energy of excitation 25(P' exceeds 
a^m^c? . This condition is equivalent to es ^ 1. Since this criterion is new to the case with 
(5 > 0, Fig. H displays the cutoff in the enhancement at > 1 for sample e^,, parameters. 

For / > it is possible for the /(/ + l)/r^ term to dominate the effective scalar potential 
everywhere, in which case the potential does not significantly modify the wavefunction and 
there is no significant Sommerfeld enhancement. This occurs for a Yukawa potential if 
+ 1) ^ l/^<^; but in the presence of a finite mass splitting this can be true even for 
60 0. 

In the €0 — )■ limit, if l{l + l)/r|. 3> e|/2, then since /(/ + l)/r^ grows faster than V{r) 
as r — )• 0, this term dominates for all r < vt- For r > r^, the effective scalar potential 
becomes y(r)-^/e|, which for a Coulomb (or Yukawa) potential falls faster than as 
r — )• oo. Thus if /(/+1) ^ (esrx)'^ /2, the potential term is always subdominant; only partial 
waves with /(/ + 1) < (esTx)'^ /2 experience a large enhancement (note that since r-p < 
this condition requires that /(/ + 1) < 2/ej). Similarly, if Vr either /(/ + l)/r^ > V{r) 
or el ^ ^(^)) the enhancement also cuts off. This is a more stringent constraint than 
the previous one when e„ > e^; if V{rv) = e^, then a significant enhancement requires 
+ 1) ^ i^vTvy (which in turn requires /(/ + 1) ^ 1/e^)- 



w 10 




0.01 0.10 1.00 10.00 



Figure 1: The enhancement S as 
a function of eg for = 0.2, = 
0.1, demonstrating the cutoff at large 
es ■ S was computed numerically using 
Mathematica. 



More generally, if 7^ 0, the constraints above 
generalize to the requirement that for the Ith partial 
wave to experience a significant Sommerfeld enhance- 
ment, the potential term must dominate the eigen- 
value for at least some range of r < rj-: if the poten- 
tial term is subdominant at r ~ r t then the effective 
scalar potential falls faster than the /(/ + term 
at larger r. 

In summary, a significant enhancement requires 

and only partial waves with /(/ + !) < !/£</), and satis- 
fying the condition V{r) > max(|e^ — e|/2|, /(/+l)/r^) 
for some r < can be enhanced. 



4. An approximate solution to the Schrodinger equation 

To obtain an approximate solution to the matrix Schrodinger equation in the form of Eq. 
|3.2| , we note that if the Yukawa potential V{r) is replaced with an exponential potential 



V{r) = Vqc"^'', then Eq. is exactly solvable in the s-wave case [37]. The relevant 
properties of this solution are described in Appendix^. For r > l/e^p, the Yukawa poten- 
tial can be approximated well by an exponential potential with an appropriate choice of 
parameters; the resulting solution for the wavefunction will be described in detail in § |4.2| . 

At r < !/€(/), the Yukawa and exponential potentials behave very differently. However, 
if the approximation of r-independent diagonalization discussed in §|3| holds over this region 
- i.e. if rx ^ l/e^ ~ then solving the matrix Schrodinger equation in this region reduces 
to solving a pair of scalar Schrodinger equations with attractive and repulsive Yukawa 



potentials. This problem will be discussed in detail in §4.1: in summary, we will show that 



the WKB approximation is valid (except possibly for an isolated turning point, denoted 
r = r*, in the repulsive case) for r > 1, breaking down when V{r) < (which occurs 
at r ~ l/e^), or when the diagonalization approximation fails at r ~ r-p. For r < 1 the 
Yukawa potential can instead be well approximated by a Coulomb potential (since we have 
assumed £</)<!) and the scalar Schrodinger equation can be solved exactly. 

We will define r = rjv/ to be the radius at which the WKB wavefunction (valid for 
1 ^ ^ iiiiii(l/e()i, r^)) is matched onto the solution obtained by approximating the Yukawa 
potential with an exponential potential (valid for r > l/e<^)- In this work we will choose 
the matching radius vm such that V{rM) = iiiax(e|/2, ep: this ensures that the WKB 
and diagonalization approximations are both valid for 1 < r < vm (except possibly at 
some isolated turning point). Fig. |2| summarizes the different regions and the approaches 
employed to obtain approximate wavefunctions in each regime. 
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Figure 2: Example of the different r-regimes and matcliing points r ~ i,r*,rM, for a sample 
parameter set (e^, = O.f, es = 0.02, — 0.05). Top panel: The ground-state (red) and excited- 
state (blue) components of the eigenstate tp- (Eq. 5^), as a function of r. When the eigenstate 



components are nearly constant, as for r < tm (< rr), the eigenstates evolve independently and 
the problem reduces to solving two decoupled scalar Schrodinger equations, r > tm contains the 
non-adiabatic transition region discussed in Middle panel: The Yukawa potential {solid black 
line) and the approximate potentials we employ, in their regimes of validity. In the r < region 
where the eigenstates are decoupled, for r < 1 the potential dominates and is well approximated by 
V{r) « 1/r [red dotted line), whereas for 1 ^ r < the WKB approximation is employed to obtain 
an approximate wavefunction. At r > vm the WKB approximation may break down, but there 
V{r) K, Voe~^^ {dashed blue line). Bottom panel: The eigenvalue A+ (red), which passes through 
zero (signalling a breakdown in the WKB approximation) at r = r*; the WKB wavefunction must 
be continued across r = r* as described in 54.1. 



The criterion that ?'t ^ 1 /^(j) is equivalent to requiring W{2e^/ej) > 1 (here W denotes 
the Lambert 1^-function) , which is true provided that e^f, > £5/2. Provided that this 
condition holds, the non-adiabatic transition described in §^ occurs close to the range of 
the interaction, where the Yukawa potential begins falling steeply. If this condition does 
not hold, then the Yukawa potential is still close to a 1/r potential when the non-adiabatic 
transition occurs, and is not well approximated by an exponential potential for all r > r^. 
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Thus our results may be less accurate in this regime, although if the kinetic energy term 
dominates for all r > vt, which is true if ^ e^, then the poor approximation of the 
true potential is irrelevant. For a Coulomb potential, a better approximate wavefunction 
might be obtained by a straightforward extension of this method: approximating the true 
potential by an exponential potential only for r ~ rj^, and matching onto a known solution 
for r >rx- However, the condition > e|/2 corresponds to requiring 5 < am^, which is 
generically satisfied for the phenomenologically motivated models of greatest interest p^ , 
and in this part of parameter space our method is very accurate. 

4.1 The small-r solution 

For V{r) ^ £5/2, as previously discussed we take the solution to the matrix Schrodinger 
equation to be (?i)+(r)^+(r)+(^_(r)^/;_(r), where ip±{r) are the energy eigenstates and (t)±{r) 
are solutions to the scalar Schrodinger equation (t>±{r) = (/(/+l)/r^— e^+e|/2ibl/(r))(/)-|-(r). 
Let us first specialize to the s-wave case with I = 0. 

We will employ the WKB approximation to find approximate wavefunctions in the 
small-r region. When the constant term — + e|/2 dominates, which is possible even if 
V{r) > e^/2 provided e^, > e^, the WKB approximation holds trivially. In the potential- 
dominated regime, \±{r) ~ ^V{r), the validity of the WKB approximation is determined 
by the parameter, 

\VWi/y{r)\ = (l/2)^A^e^W2(e^ + 1/^). 

For r <^ l/e^, the approximation is good so long as l/2y/r < 1, and breaks down at r ~ 1. 
At r ~ where the potential begins falling rapidly, the condition for validity of the 

WKB approximation becomes e,f,/\/V{r) < 1, so the approximation breaks down when 
V{r) ~ e^. By the assumption of a significant enhancement, l/e^ > 1, so there is an inter- 
mediate range of r where the WKB approximation is valid. At r < 1 we can approximate 
V{r) ~ 1/r and assume the potential term dominates (since ^ 1 by assumption), solve 
the resulting ODE exactly, and match onto the WKB solution for r > 1. 

Note that we can now see explicitly that if > 1 the enhancement must cut off, since 
this corresponds to the potential falling off sharply at r < 1, and in the r ~ — 1 region 
the magnitude of the wavefunction does not evolve significantly. Similarly, if e^, > 1 then 
the kinetic term begins to dominate at r < 1, and the wavefunction is well approximated 
by a plane wave from that radius outward. If > 1 then r-r ^ 1, and at r > the 
effective attractive scalar potential for the ground state is given by ~ l/r^ej (Eq. |3.5D (for 
the excited state the sign of the effective potential is reversed). While the potential term 
dominates (and r > vt), the approximate solutions to the Schrodinger equation scale as 
^i/2±-y/i/4-i/e^ (^Qj. j,i/2±-^i/4+i/e| £qj, excited state). The ground state wavefunction is 
therefore oscillatory provided es < 2, with the amplitude scaling as l/\/r, but for es > 2, 
the physical solution becomes nearly constant (exactly constant in the limit of large e^). 
Thus the amplitude of the wavefunction does not grow significantly for r > in this case, 
and as stated previously, no large enhancement is possible for » 1. 

Returning to our calculation of the small-r wavefunctions, the eigenvalue A+ is negative 
when the kinetic term dominates but positive when the potential term dominates; if ey > 
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es, it may possess a zero in this small-r region {V{r) > £5/2)- In this case the WKB 
approximation for the "repulsed" solution has a turning point in the inner region, which 
we must take into account. 

Let us first consider the potential-dominated region r < 1, where the WKB approx- 
imation breaks down and we can set X± ~ itl/r (since by assumption e^,, < 1). The 
solutions to the scalar Schrodinger equations are then given simply by, 

(p-ir) = A_y/^Ji{2y/7)-7r^_{0)y/7Yii2y/^), 0+(r) = ^+\/f/i(2V^)+2(/.+ (0) V^i^i(2V^). 

(4.1) 

At r > 1, using the large-r asymptotics of the Bessel functions, these solutions match 
smoothly onto the WKB solutions. 



(4.3) 



Let r* be the radius such that V{r*) = e^y e'^ — ej, so A+(r*) = (in the case > es). 
Then linearizing the potential across the turning point and matching to the WKB solutions 
on either side, we obtain the WKB solution for r > r*: 



We can check this result in the limit — ?■ 0, — ?■ 0: we should recover the usual 
expression for the Gamow-Sommerfeld suppression factor due to a repulsive Coulomb po- 
tential n, m, 

s ^' ' 



In this case r* = and f^* y^J^dr' = 7r/2e^. Imposing the boundary conditions that 

(j)+{0) = 1 and (j){r) is purely outgoing at r — )• 00, we find, 

^+ = ^^/ W--- W|<^+(oo)| 



2/ " Ve„e-A--i' 



(note that only the integral of ^y X-^ from to r* contributes to the amplitude, as beyond 
r* the integrand is imaginary), which yields a suppression factor, 



ey Ye'^A" - 1 + (l/4)e-'^A" 



Dropping the exponentially small term in the denominator (since we have assumed that 
<C 1 in the matching with the r < 1 solution), we obtain the correct expression. 
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For the higher partial waves which experience a significant enhancement, there must 
exist a region where V{r) ~ 1/r and is large compared to both the l{l + l)/r^ term and 
the term. For / = 0, this region extends inward to r = 0; for higher partial waves the 
/(/ + l)/r^ term will dominate for r < Z(Z + 1). In the range of r where V{r) ~ 1/r and 
either V{r) or /(/ + l)/r'^ is the dominant term, we can approximate the eigenvalues as 
A-t = l{l + l)/r^ lb V{r), and solve the Schrodinger equation exactly. For larger r, but still 
in the small-r regime r < l/e^, the l{l + l)/r^ term is subdominant (since it falls faster 
than V{r)), and so we can ignore the effects of non-zero I, and match onto the s-wave WKB 
wavefunction. The discussion of the validity of the WKB approximation in the s-wave case 
is equally applicable to this case: the WKB approximation may break down for r < 1(1 + 1) 
where the l{l + l)/r'^ term dominates, but by the assumption of a significant enhancement, 
> /(/ + 1), so again there is a range of r where the WKB approximation is valid. 

Thus the only required change to the s-wave calculation is to approximate the eigen- 



values at low r by X±{r) = l[l + l)/r it 1/r, replacing Eq. 4.1 with 



0_(r) = r(2(/ + l)) ^_^j^^^^(2^) _ ^^^^^^ ^^ <^_(0)V^y,,+i(2Vr), 

Mr) = r(2(^ + l)) ^^^^^^^^(2^) ^ __?__<^^(0)^K2/+l(2^^^). (4.5) 

Here (/>+(0), i^-(O) are the coefficients of the r~' terms in (j)^{r), 4>~{r) respectively, as 
r — 0. The large-r limits of these solutions match onto the WKB wavefunctions derived 
for the s-wave case, with the replacements (/'±(0) ci)±{0)/T{2l + 1), A± A±r{2l + 2)/2. 
The WKB wavefunctions also gain an overall multiplicative factor of (—1)' (note that the 
Q±Jv^dr i^gj-jjig in the WKB wavefunctions should still be taken to refer to the s-wave 
eigenvalues X±). 

4.2 The large-r solution 

Where the WKB approximation breaks down (at V{r) ~ e^) and/or the diagonalization ap- 
proximation fails (at V{r) ~ e^/2), we match the previously derived wavefunctions onto the 
exact wavefunction in the presence of an off-diagonal exponential potential V{r) = VQe~^^ , 
outlined in Appendix ^ For convenience, we define a new parameter z = V^j^e^^'^^/lG/i^. 

The parameters of the exponential potential, Vq and fi, are set by the requirement that 
the exponential potential should mimic the Yukawa potential for r > vm, the matching 
radius. We follow [^] and require that rV{r)dr = rV[r)dr We also impose the 
condition that V{rM) = ^(^m)- 

The parameter is then given by. 



- = -{rM + -\^ ^i = eA\ + \Jl + ^\. (4.6) 
60 /^V ^ \2 2\i e^TM / 



^This criterion comes from solving the Lippman-Schwinger form of the Schrodinger equation, to obtain 
Ri{r) as a series in the potential V\ in the limit e*, e„ — > 0, the condition that Ri{0) is correct to 1st order 



in the coupling a fixes J rVdr 
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A more sophisticated prescription for the potential matching might yield better agreement; 
however, we will show that this simple prescription already gives good agreement with 
numerical results. 



We will employ Eq. 2.16 to determine the Sommerfeld enhancement; thus, we are 
interested in wavefunctions which are purely outgoing (or exponentially decaying) as r — )• 
oo. As discussed in Appendix |A|, this condition implies Ci = C3 = in Eq. A.2| . The 
remaining coefficients C2, C4 are determined by matching to the solutions in the inner 
region. For ease in reading off the Sommerfeld enhancement, we write the exact solution 
for the exponential potential in the form, 
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Here the p-Fg's are hypergeometric functions. The elements of the Sommerfeld enhancement 
matrix T relevant to particles initially in the ground state are then given by rj (calculated 
for the two LI boundary conditions at r = 0, ip{0) = (0,1) and il^{0) = (1,0)), and for 
particles initially in the excited state by (. For et, < e^, only the r] coefficient is relevant, 
as the two-body excited state cannot be on-shell. 

The matching procedure in the inner part of the transition region is described in 
Appendix p. We define: 



X± = -el + ej/2±J{ej/2y + V^, 



r^: the radius (if any) at which the eigenvalue for the repulsed eigenstate of the 
exponential potential passes through zero, defined by V{r* 



<5' 



rs- a (possibly negative) quantity chosen such that Vqc ^> e^,e^. The corre- 
sponding value of z is denoted zs = z{rs)- 



• CA = 

We solve for ij, C for arbitrary boundary conditions at the origin, obtaining for s-wave 
annihilation: 
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Below threshold: 



-(0) 



Above threshold: 



-1 + e" 



1 _ 


r 







2/x 2^t 2 



EA _|_ 1 

" 2^j 2/^ 2 



e 



1 _ 


r 







i(j)-{0)e M 



2m y 



r 

(O)e' 



2/x 



+ 



-1 + e 



-+2je- 



26** 



v 



-1 + e" 



-2^/2-2 



1 _ HA 



2/^ 



X li<^_(0)e^+^^- + 



r 

(0)e 



£A— 
2ft 



g A. _|_ g M 



-1 + e 



1 — e A" 



2 - g2(eT-0y) + 2e2(^*+^T-ey) + g-2e* 
In both the above and below threshold cases, the phases are given by, 



(4.7) 



/ JX-dr-Aizy- ^XZdr, 6+= Jx+dr-Az^- 

Jrs Jo Jrs Jo 



rM 



X^dr, 
(4.8) 



&T = J yX+dr, = J vA+dr, = J ^JX+dr. 

If Z > 0, then there are two cases of interest. First, suppose that the /(/ + l)/r'^ term 
is subdominant for all r > tm- Since the potential term is rapidly falling for r > r^f, this 
requires that /(/ + l)/r\i < el- Then the large-r s-wave wavefunction described above can 



be matched onto the rescaled WKB solution outlined in §4.1, and the only change from 
the s-wave case is that the elements of the enhancement matrix T must all be rescaled by 
a factor (— l)'/r(2/ + 1). The Sommerfeld enhancement for / > is then identical to the 

s-wave case, up to an overall rescaling by ^ r(2/+i')e' ) (^sing Eq. |2.16| ). 

If e,„ < + l)/rM, this argument no longer holds, as the /(/ + l)/r'^ term cannot 
be neglected for r > tm- Since the main focus of our work is the s-wave enhancement, we 
will not treat this case here. 



5. Properties of the Sommerfeld enhancement in an inelastic model 

The results derived above can be used to compute the Sommerfeld enhancement for any 
annihilation matrix T, and particles in any initial state. However, in most applications 
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we will be interested only in particles initially in the ground state. Furthermore, in the 
simplest models, all elements of the annihilation matrix F are identical. 

In this case, the s-wave annihilation enhancement for particles initially in the ground 
state (Eq. |2.16 ) reduces to, 

S =\Tn+Ti2\^ = \Moo)\^ with V(0) = (1,1). (5.1) 

This case is therefore the closest 2-state analogue of the purely attractive Yukawa potential, 
as the boundary condition at the origin sets the repulsive component (/>+(0) to zero (and 
(j)-{0) = \/2). The expressions derived in the previous section simplify dramatically in this 
case, yielding, 



S = — sinh 



(e„7r//i) -COS (^-y/e2 -e2^/^+26»_) 



cosh 



cosh ( (e„ + ^-e2+e2^ n/2fi'^ scch ( (e„ - ^-ej+e^'^ 7r/2/x) 
cosh((e„ + ^-e2+e2^7r//i)-cos(2e_) 



(5.2) 



While this expression still contains a numerical integral for 6^ , it is well-behaved and 
fast to compute ^. 

More generally, note that whenever 0+(O) appears in Eq. 4.7, it is accompanied by 
a factor of e^+. When < es, 0+ is purely real (since A+, A+ are both positive for any 
r G [0, ta,/]), and can be written in the form, 



J mi 



^ ^ 2 



-2\ 2 



+ le/i^z - 2[xz^!^ dz 



'max(16e^,4e^) 



Here zs has been taken to oo. The integrand in the first term satisfies the inequality, 

= \ 



1 



2^z 



\ 



ei + -f + 



2\ 2 



+ IQix^z - 2nz 



1/4 



(5.4) 



Here we have used the relations \/a + b < y/a + Vb and \/l + a < 1 + a/2 foi a, b real and 
positive. 



Note that the J^'^'' -J |A_ \dr — 'iz]J'^ contribution to 9- may be best computed by changing variables in 



the integral to z, and rewriting the second term as an integral, giving 
1 



I \ " 2 



-2\ 2 



+ l&ii-^z - 2nz 



1/4 



dz — 



/max(e^,4/4)/16M4 

This prevents the integral from diverging at large 2, so 25 can be taken to 00 



max(16e0, Aej] 



1/4 



Since the last term in Eq. ^.31 is strictly negative (the integral is always real and 
positive), integrating Eq. 5A yields, 

0^ < (^^j (4(max(4,4/4)/16/)-V4) _ f ^ < 0. (5.5) 

So in the case where e„ < e^, 9+ is real and negative. When ey > eg, the real part of 
64- becomes, 



s 

rs 



In this case we can write. 



V^=^-eS + | + y(|) + V{rf < ^-el + e] + V{r) < 



and so. 



Jts ^ Jrs ^ ^ 

Thus the real part of ^4. is always negative, and the contribution of the repulsed 
component is exponentially suppressed. To a first approximation, then, the elements of 
the enhancement matrix T are determined solely by the coefficient of (/)_(0). Since both 
the ■0(0) = (0,1) and tp{0) = (1,0) boundary conditions correspond to (/>_(0) = l/\/2, the 
s-wave enhancement in this case (for particles initially in the ground state) becomes, 

s = 5att(rii + ri2 + Tai + r22)/4rii. 



where Satt is given by Eq. 5.2 



5.1 Numerical tests of the approximation 

In Fig. ^, we compare S computed using Eq. 5.2 to the exact numerical solutions for 
the enhancement (computed using Mathematica), for = 0, 10^^ 10^^, in the regions of 
parameter space where Eq. 5.2 is expected to hold and the numerical solution is stable. In 
all cases the approximate form correctly reproduces the numeric results to good accuracy, 
although it becomes less reliable as e^, e„ — )• 1, as expected. 

5.2 Magnitude of the non-resonant enhancement 

For e„ ^> e^, e^, ^> fi, we obtain, 



evj cosh -cos(26'_) 



and furthermore the exponential terms dominate the oscillatory cos term, yielding S ~ 
vr/e^. Thus we recover the usual Coulomb enhancement in the limit where the mass splitting 
and mediator mass parameters are small compared to the velocity parameter e^,, as required. 
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Figure 3: The Sommerfeld enhancement as a function of and e„, computed both numerically 
and using the semi-analytic results derived in this work, for (top) es = 0, (middle) eg ~ 0.01, 
and (bottom) es = 0.1. For purposes of comparison, we restrict the range of to regions where 
the numerical solutions to the Schrodinger equation are stable. Left panel: The enhancement 
computed by numerically solving the Schrodinger equation using Mathematica, Middle panel: The 
enhancement computed from Eq. ^.2| , Right panel: the S — 10, 100, 1000 contours for the numeric 
calculation of S (solid black line) overlaid on the same contours for the approximate calculation 
(dashed red line). 



If et, 3> /u but < e^, on the other hand, then as previously the exponential terms 
dominate the oscillatory term, but now ~ 27r/e^. So below the threshold for on-shell 
excitation of the upper state, the non-resonant enhancement is increased by a factor of 2. 
This effect is illustrated in Fig. |^. 

When ey < fi, then up to a factor of the form 1/(1 — cos 6*) describing the resonances 
(which will be discussed next), both above and below the excitation threshold, the enhance- 
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Figure 4: The Sommerfeld enhancement as a function of e^, for — 0.01, and {solid line) 
es — 0.1, (dashed line) es = 0, using the approximate expressions described in the text. Left 
'panel: The left-hand side of the plot (e^ ^ Ct, < ea) demonstrates the increase in the non-resonant 
enhancement discussed in § ^.2| , whereas the right-hand side [t^ < e^) shows the higher, shifted 
resonances discussed in § ^.3[ Right panel: As the left panel, but also displaying the results of 
numerically solving the Schrodinger equation, in the region of numerical stability. Black diamonds: 
es = 0.1, black dots: es = 0. The approximate solutions reproduce the numerical results well, and 
the predicted larger resonances and shifts in the resonance positions are clearly visible. 



ment saturates at the same value S ~ 2tt / fi. In the regime where Eq. 5^ is expected 
to be accurate (r^ ~ l/cc^)) A* ~ e</>5 and this is the usual saturation of the enhancement 
occurring when drops below e<^. As e,^ — )• 0, /x — )• (this holds regardless of es, provided 
only that < 1), and so the enhancement does not saturate (scaling as l/e^, down to 
arbitrarily small velocities) ^. 

5.3 Resonance behavior 



When falls below e^, the resonance structure described by Eq. p.2| changes in two major 
ways. First, the heights of the resonances change: where < e^, the resonances can be 
significantly larger than in the zevo-S case. This can be seen simply by examining the 
cosh term in the denominator of Eq. |5.2| , where <^ fi: the height of each resonance at 
its peak is controlled by l/(cosh(et,7r/;u) — 1) « (/x/vret,)^ in the below-threshold case, and 

l/(cosh((e^-|- Y^e^ — e^)7r///) — 1) f» (/i/7r(e„-|- y^'e? — ^5))^ above threshold. Thus we expect 
the resonances to be ^ 4x higher in the below-threshold case compared to the zero-5 limit, 
as illustrated in Fig. ^ 

The resonances also shift to different values of when e^ falls below eg. We can read off 
the resonance locations from the semi-analytic solution: above threshold they occur when 



*While our method may be expected to become inaccurate for < eg/2, this general statement - that 
in the = limit, at low velocity S scales as l/e^, independent of eg provided < 1 - is expected to 
remain valid, based on matching the WKB solution in the small-r region onto the exact solution for an 
effective 1/r^ potential in the large-r region. 
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9^ = nvr, n £ Z, whereas below threshold they shift to satisfy 6- — ^rye^ — = nvr, 

n G Z. To estimate the resonance locations we need an analytic estimate for the integral 
0-. 

In the range 1 > > e^, et, < /x, 6*- is very well approximated^ by —sJI-Kje^. The 
condition that Ct, < /i is required for significant resonances (since otherwise the oscillatory 



terms in Eq. |5.2| are entirely dominated by the exponential terms), and since the height of 
the resonances is proportional to /i^ ~ e^, the most pronounced resonances occur at large 

where this simple estimate for is valid (for < there are significant corrections 
to this expression, but no especially simple form has been found). 

Employing this estimate for we find that above threshold, resonances occur where 

~ {2/t:)/v? , n G J\f. Below threshold, we write /x = ce^ where c is a slowly varying 0{1) 
function of e^, and approximate c = constant. Solving the resulting quadratic equation for 
the resonances at = 0, we obtain. 
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1 [2 nes 2 n j-\ 

2n^ \7r c vr / 



It appears that this equation imposes an upper bound on n: however, recall that 
we have assumed e,^ > e^. If the term under the square root is to approach zero, this 
corresponds to n ~ c/vre^, and ~ 1/(27™^^^^^^) ~ (7r/2)e|/c^ < e^. Thus by virtue of our 
previous assumptions, we cannot probe the large-n limit: we will restrict ourselves to the 
case where nvre^/c ^ 1. In this regime, we can approximate the solutions for the resonance 
positions as: 

2 1 es I (,____ ,_,2^ 



O {mres/c) 



Only the first set of solutions lies in the > es regime that we have assumed, so the 
effect is to shift the existing resonances downward in e^f, by es/cn (note that relative to 
the resonance position, this corresponds to a fractional shift of Tmes/2c, i.e. the fractional 
shift grows with increasing n). Fig. Q clearly shows the increasing fractional shift to lower 
values of e^f, with increasing n. 

5.4 Behavior at the = es threshold 

In the 5 = limit, is a monotonic function of e^,. This no longer holds true when (5 7^ 0: 
the position of the resonances varies as a function of e^, and this effect can give rise to 
noticeable spikes in the enhancement around the threshold velocity = (where the 
resonance positions are most dependent on for far from 65, the resonance positions 
are largely independent of e^). In particular, the Sommerfeld enhancement near threshold 
may be even larger than the saturated enhancement at low velocity. In the context of DM 
annihilation, this effect could alleviate tension with early-universe constraints on the rate 
of dark matter annihilation (e.g. 41, 42, Q), since in the early universe the DM 



is extremely slow-moving, while in the present-day Galactic halo is of the same order as 



^This expression was derived by noting that in this regime the j^^' y^X-dr contribution to 0- dominates, 
and setting A_ « —V{r); the estimate was then checked numerically against the full integral form of 9-. 
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€5 for phenomenologically interesting mass splittings ®. To investigate this possibility we 
examine the ratio, 



_ S{ey = es) _ fi 

-'•'thres 



sinh (2^) (1 - cos + 20.{e, = 0))) 



S{ey = 0) Ties cosh j - cos (20_ (e„ = e^)) 

If ires > /i, then the exponential terms dominate, and i?thres ~ (t^/'^^s) ^1 — cos + 26 -{e^ = 0)^^ • 
We see that i?thres < 2^/ ires < 2, and the maximal value of -Rthres will occur when /i ~ Tre^ 
and ires I n + 2^_(e^ = 0) is as close as possible to (2n — l)7r, n G M . An example where 
i?thres attains a large value is shown in Fig. |5[ Of course, in realistic applications it will be 
necessary to integrate the enhancement over some velocity distribution, which will tend to 
smooth out the resonance: the effect of smoothing by a Maxwell-Boltzmann distribution 
is shown in the right-hand panel of Fig. ^. 




Figure 5: Examples of non-monotonicity of the Sommerfeld enhancement as a function of c^. 
Left panel: the enhancement computed by numerically solving the Schrodinger equation, for (solid 
line) es = 0.01, = 0.04, {dashed line) eg = 0, = 0.043. The slightly different values of 
£0 are chosen to give the same behavior in the et, — > limit, to show the qualitatively different 
behavior at intermediate velocities for models which are similar at very low and very high velocities. 
Middle panel: the enhancement computed using the approximate results derived in this work, for 
{solid line) eg = 0.01, = 0.0414, {dashed line) eg = 0, = 0.0433. The marginally different 
parameter choices for the analytic approximation and the numerical calculation reflect the small 
degree to which the approximation shifts the resonances; in both cases the parameters were chosen 
to demonstrate a clear enhancement at e„ ~ e^. Right panel: the effect of smoothing the velocity- 
dependent resonance (using the numerical results, and the same parameters as the left-hand panel) 
with a Maxwell-Boltzmann distribution e^e"*^"/^*^" . The enhancement is plotted as a function of 
the dimensionless velocity dispersion e^r- 



In the opposite limit, where vre^ <C /i, 



l-cos(^^ + 20_(e^ = 0) 

-^thres ~ ~7 2 

-cos{2e.{e, = es)) 



®For example, for a 1 TeV WIMP with 100 keV mass splitting, e„ ~ es corresponds to w/c ~ 130 km/s, 
which is very close to estimates for the one-particle rms velocity of the WIMPs |45[. 
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There is the potential for a large enhancement if cos(20_(e^ = e^)) ~ 1, leaving only a 
small term in the denominator. However, at lowest order 9- is independent of e^,, so the 
numerator would then be zero at lowest order. Consequently, in order to determine i?thres 
in this limit we must compute the next-order term in 0_ for e^, = and = e<5. 

Away from the transition region (where V and V cancel out in ^_ ) , we can approximate 
y/X^ = iVV ±iej/4:VV, where the + sign corresponds to = es and the — sign to e^ = 0. 

An identical expression holds for ^/X^ with the replacement V ^ V. Thus, 



e 



5 



9_{e, = es) - 9-{e, = 0) = ^ i^J l/VVdr-J^ l/VVdr 

The first integral evaluates to 2\/^/{^es) (since F(rj') = e^/2), whereas the second scales 
like e . ; assuming n ^ e^, the first term dominates provided es < ^/^■ Thus we have. 



and consequently, 



Rthr 



.{ey = es) « 9^{e^ = 0) + V2es/fi, 



-cos( ^"-y)-^ +2g_(6. = 6^) 
i+(^Y -cos{29.{e, = es)) 



We see that as 0, i?thres — ^ 1, and if cos(20_(e^ = e^)) — )• 1, i?thrcs (1 — 2\/2/7r)^ < 1. 
As well as these limits, it can be numerically checked that everywhere in this region of 
parameter space -Rthres < 1- Consequently, -Rthres is robustly bounded above by iZthres < 2. 

5.5 The £5 —7- limit and comparison to previous results 



In the case where the mass splitting is zero, = 0, Eq. 5.2 simplifies to 



S= — 



2n\ «inh ^ cosh /^x sinh^2.e. 



/ cosh ( ^ ) - cos(26l_) \^^J cosh ) - cos(26l_) 



The phase 9- is approximately given by 0_ « —y^27r/e^ in the — )• limit. 

An approximate solution for the Sommerfeld enhancement due to a Yukawa poten- 
tial was derived in |34], using the exact solution to the Schrodinger equation with the 



Hulthen potential |T9|. In terms of the dimensionless parameters used in this work, the 
expression for the s-wave enhancement derived in can be written as. 



TT 



sinh 



S={-] ; ; •• , . -■ (5.6) 

V^»/„,h(^)-c„=(2.y^ 

where is a parameter of the exactly solvable potential used to approximate the Yukawa 
potential, and set to be = (7r^/6)e^ in |34|. We see that the two expressions are identical. 



with the replacements o and 0_ -H- vr. — At low velocities the resonance 
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positions are set by the argument of the cos term, which is ~ 27r/ y^vr^e^/G = 2\f^l ^Te^ 



in Eq. 5.6 and ~ I^JIt: j ^Jt^ in our approximate solution, in close agreement. Taking 
^'Af ~ in Eq. yields [i ~ (1/2)(1 + -v/5)e0, which is close to the 7r^e(^/6 value employed 
in (however, ///e,^ is only approximately constant, and as varies some disagreement 



between our formula and that of will ensue). The solution of [34| did not make the 
approximation that e^^e^ <^ 1, and thus is superior in the region where these assumptions 
break down and the enhancement is small, but the two approximations are very similar 
(and both reproduce the numeric results accurately) whenever the enhancement is large. 

6. Conclusion 

We have developed a detailed approximate expression for the s-wave Sommerfeld enhance- 
ment in the presence of a mass splitting, in the case where only inelastic scattering between 
the interacting particles is possible at tree level. The 5 — )• limit of our results yields an 
accurate approximation for the enhancement due to an attractive scalar Yukawa potential. 
The enhancement cuts off at low velocity if the mass splitting 6 > a'^'m^^, the effective Bohr 
energy of the two-body state, but even for smaller mass splittings, interesting modifications 
of the enhancement relative to the 5^0 case can occur. Additionally, we have derived 
simple expressions to facilitate computation of the Sommerfeld enhancement due to an 
arbitrary N x N matrix potential, where the Ith. partial wave dominates. 

In the special case where all elements of the annihilation matrix are identical, our 
expression for the s-wave enhancement takes a simple form, 

1 

cosh(€^7r//^)— cos^-y/ e^—e^Tr/ 11+26-^ 



5=^sinh(^ 



cosh ( (e„ + 7r/2Ai) scch ( (e„ - ^-ej+e^^ 7r/2Ai) 

cosh((e„ + -y/-e2+e2)7r//i)-cos(20_) 

where fi and 0_ are real functions which we have defined. This expression is valid pro- 
vided rriff) < am-^, v/c < a, and S < Q^m^, although it may become less accurate for 
6 > max(am(^, m^{v/cf'). If any of the three former conditions are violated, there is no 
significant enhancement {S ~ 0{1)). In the regions where this expression is expected to 
be valid, comparison to numerical results demonstrates that it accurately reproduces all 
significant features of the enhancement. Furthermore, while derived for the case where all 
elements of the annihilation matrix are identical, a simple rescaling of this result by the 
average of the elements of the annihilation matrix (normalized to the Fn element) can 
accurately describe the enhancement for a general annihilation matrix. 

Employing this approximate form for the enhancement, we have shown that below the 
threshold for on-shell production of the excited state, in comparison to the case with no 
mass splitting: 

• The non-resonant, unsaturated enhancement is higher by a factor of ~ 2. 

• The positions of the resonances shift to lower values of {m^/am^) and become 
more widely spaced. The heights of the resonances increase by a factor of ~ 4. 
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• The enhancement S is no longer a monotonic function of velocity: close to the thresh- 
old for on-shell excitation of the higher-mass state, the enhancement may be as much 
as a factor of ~ 2 greater than its saturated value (i.e. its value in the limit u — t- 0). 

In the context of DM annihilation, if the kinetic energy of WIMPs in the neighborhood 
of the Earth is close to the threshold for excitation of the higher-mass state (which is true 
by construction for models which realize the iDM mechanism), increased near-threshold 
enhancement may slightly alleviate constraints on DM annihilation from the early universe, 
where the DM is extremely slow-moving. The generically larger values of the unsaturated 
non-resonant enhancement, and the shift in position of the resonances, may help provide 
sufficiently large enhancements to generate recently observed cosmic-ray excesses, particu- 
larly for lower-mass force carriers. 
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A. The exact solution for an exponential potential 



We will employ an exactly solvable exponential potential to analyze the non-adiabatic 



transition at V{r) ~ e|/2. If we take V{r) = Vqc then the Schrodinger equation 



i^iir)] ^ / -el -V{r)\ Mr) 
i^'iir) \-y{r) ej-el {Mr) 



can be solved exactly by the substitution |37| z = VqC ^'^''/16^^, since then the 4th-order 
ODE for takes the form. 



' ' d 



n=l 



Yl(z—-bn)i'l + zi,i= 0, (A.l) 



6i = iej2^, 62 = -iev/2fi, 63 = 1/2 + i^e^ _ ^2/2^, ft^ = 1/2 _ i^e^ - 

which is exactly solvable in terms of hyper geometric functions. V2 can be deduced imme- 
diately from Tpi, since ip2 = {ipi + e?V'i)/^(^); a^id simplified using the hypergeometric 
function identities, 

(z-^+bk-l) pFg{ai, . . . ,ap;bi, . . . ,bk, . . .,bq;z) = (bk-l) pFg{ai, . . .,ap]bi, . . . ,6^-1,. . .,bg;z) 
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^ pFqiai, . . . , Op; 61, . . . , bq; z) = ^| pFqiai + 1, . . . ,ap + l;bi + 1, . . . ,bq + 1; z). 

We find, 



V'l = (^22 ^'^ 0^3 



<=2 -U f2 
^ 



2ji ' 2 2/i 



2/x 



^^'"^^ m'2 2/x 2/x '2^2//^ 



2// 



+ VZ C4Z 2m 0F3 



^ 2 2n 2n ' 2 2fi 
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■^5 ^ ""v 
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■^5 ^ ''v 



1p2 = -^fz 



C2Z 2m 



(2 ^2^ 



2 .2_,2 
+ 



2ix) ^ in 



ClZ2M 



'1 _|_ -L 

^2^*2^^ ^ 



0-^3 



0-P3 



3 ie„ + 3 ie,; 



2 2;u 2iJL ' 2 2n 



+ 



^5 + 



2/x 



2 2/x 



2(1 



"' 2 + 2/x + 



c2 _L ^2 



2// 



-CiZ 2m 



X oi^3 



{}, l- 
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/i ' 2 2/x 



u 1 ^ ze,; 
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^5 



1 



C^Z 2M _1L + + i 



2/x 



2^ ' 2 2fi 



J 



2n 



X 0-P3 



1 + 



-2 _ ^2 



"'2 2/x ^ 



-U f2 , 



2^ 



- - A ^ + 

' 2 ^ 2/i ^ 



2ju 



,1 + ^ 



,2; 



(A.2) 



The hyper geometric functions asymptote to 1 as r —t- 00, so the prefactors determine 
their asymptotic behavior. We see that the four linearly independent solutions with coeffi- 
cients Ci, C2, C3 and C4 correspond to outgoing and ingoing spherical waves in the ground 
and excited states as r ^ 00: 



Ci term: ingoing wave in the ground state, suppressed by ^ at large r in the excited 
state. 
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C2 term: outgoing wave in the ground state, suppressed by ^/z at large r in the 
excited state, 



• C3 term: ingoing wave in the excited state, suppressed by ^/z at large r in the ground 
state, 

• C4 term: outgoing wave in the excited state, suppressed by ^/z at large r in the 
ground state. 

This solution can immediately be employed as an exactly solvable toy model for the 
Sommerfeld enhancement (similar to the finite square well, as mentioned in but the 

resulting expression for the enhancement does not have a particularly enlightening form, 
so we do not discuss it further in this work. Instead, we employ this exact solution to 
describe the behavior of the wavefunction at large r. 



B. Matching of wavefunctions at r ~ tm 

In the limit z — )• 00, the exact solution for the exponential potential (hereafter termed 
the "exponential solution" for brevity) takes a relatively simple form, determined by the 
asymptotics of the hypergeometric function: 



P rn J , , 1 _ T{a)T{b)T{c) ^ 
0F3 [{}, {a, 6, c],z\^ ———j—z-i 



(B.l) 

where 7 = (1/4) (3/2 — a — 6 — c) [0. There is also an exponentially suppressed term, of 
the form cos (4^7)e-^^ ' HI; provided < ^, the imaginary part of 7 is small, and this 
term can be neglected (the effect of including this term has been tested, and it does not 
significantly modify our results). Where > the enhancement to annihilation takes a 
simple form, independent of the details of the matching procedure (i.e. terms involving 
the parameter 9- are subdominant, and the dependence of 5 on /x cancels out), and again 
this term can be neglected ^. 

To determine the form of the exponential solution near the transition region, we use 
the WKB approximation to propagate the known asymptotic (large z) solution in to the 
transition region. We can then match onto the WKB solution for the inner region (§[4.1j). 

For large z, the potential term dominates and the eigenstates are approximately 
(— l,l)/-v/2 (for the repulsed component) and (l,l)/-v/2 (for the attracted component). 
Writing the large-z form of the exact solution (Eq. |A.2| ) in the form + we can 

extract the solutions (j)± to the corresponding scalar Schrodinger equations. 

Let us write the WKB solutions in the form. 



|A±|V4 



^Note that neglecting this term may give incorrect results if these wavefunctions are employed to study 
scattering in a similar system, in the > A* regime, since in the scattering problem the phase shifts induced 
by the potential are physically significant and depend on the matching procedure. 
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where A± = -el + ± y (e^/2)2 + V{zY, and rs is chosen so that ^(rs') > e^, e| (we 
choose this form to facihtate matching both onto the previously derived WKB solutions 
for the Yukawa potential, and onto the asymptotic forms of the exact solution for the 
exponential potential) . 



Choosing some r > rg that still lies within the potential-dominated region, the WKB 
solutions match smoothly onto the exact solution for the exponential potential, yielding 
the relations: 



(27r)3/2 



(B.2) 
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(B.4) 
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As discussed above, we neglect the exponentially suppressed Ej^ term. 

Below threshold, there is no turning point in the WKB approximation for the repulsed 
eigenstate. At some matching point tm, chosen so that the WKB approximation holds for 
both potentials and ^(rjvf) ~ ^iju)-, the matching equations become, 



(B.5) 



1 



^"/^ {A. + i7r<^_ (0) ) , /3_ = e^"/^ - (0) ) , 



20F 



20F 



iA/vr 2^/7r 
Above threshold, suppose there is a turning point in the inner region, at r = r^, for 
the repulsive eigenstate of the exponential potential. Then the matching equations for the 
repulsive eigenstate are modified to give. 



1^ 



(B.6) 



1 [ iA+ 
2i \ 2y/^ 



(B.7) 



2^ 



A, 



Jo"[V^+'V^)dr-2j^- ^X+dr _ ^-2 /J'* ^dr+ J^,'^' [^-^ X+)dr 



In either case, we now have six equations (Eqs 



|B.5| and Ej^ = 0, combined 



with Eq. p.5| and/or [B.?] and |B.8| , plus the boundary conditions which fix two of the 
C coefficients) in six unknowns (Ci, C2, C3, C4, and ^-). Thus we can solve for the 
large-r behavior of the exponential solution in terms of the boundary conditions at the 
origin; the results are summarized in §4.2. 
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